

%% Load WRF output ***Cumulative variable***
%% Take precipitation as an example

tic

% Loc
Loc = '.\';

% sec
Sce = { '1995', '2018' };

% years
time_year = 2018 : 2020;

% Load
W_Pre_bef = zeros( [size(W_lon), 365, length(time_year) ] );	% UNEP
W_Pre_aft = zeros( [size(W_lon), 365, length(time_year) ] );	% CTRL

for y = 1:size(time_year,2)
   
    for m = 1 : 12

        md = eomday(2018, m);

        for d = 1:md

            n = sum( eomday(2018, [1:m-1]) ) + d;

            % bef.
            flag_Sce = 1;
            file = [ Loc, 'a_', num2str( time_year(y) ), '_', Sce{ flag_Sce }, '\wrfout_daily_d02_', ...
                    datestr( datenum(time_year(y), m, d), 'yyyy-mm-dd' ) , '_00_00_00.nc' ];
            file_tom = [ Loc, 'a_', num2str( time_year(y) ), '_', Sce{ flag_Sce }, '\wrfout_daily_d02_', ...
                    datestr( datenum(time_year(y), m, d)+1, 'yyyy-mm-dd' ) , '_00_00_00.nc' ];
            RAINNC_tod = rot90(  double(   ncread(file,'RAINNC')   )  ) + rot90(  double(   ncread(file,'I_RAINNC')   )  )*50;	% 50 is the bucket value in WRF, please be consistent with namelist
            RAINNC_tom = rot90(  double(   ncread(file_tom,'RAINNC')   )  ) + rot90(  double(   ncread(file_tom,'I_RAINNC')   )  )*50;	% 50 is the bucket value in WRF, please be consistent with namelist
            RAINNC_bef = RAINNC_tom - RAINNC_tod;

            % aft.
            flag_Sce = 2;
            file = [ Loc, 'a_', num2str( time_year(y) ), '_', Sce{ flag_Sce }, '\wrfout_daily_d02_', ...
                    datestr( datenum(time_year(y), m, d), 'yyyy-mm-dd' ) , '_00_00_00.nc' ];
            file_tom = [ Loc, 'a_', num2str( time_year(y) ), '_', Sce{ flag_Sce }, '\wrfout_daily_d02_', ...
                    datestr( datenum(time_year(y), m, d)+1, 'yyyy-mm-dd' ) , '_00_00_00.nc' ];
            RAINNC_tod = rot90(  double(   ncread(file,'RAINNC')   )  ) + rot90(  double(   ncread(file,'I_RAINNC')   )  )*50;	% 50 is the bucket value in WRF, please be consistent with namelist
            RAINNC_tom = rot90(  double(   ncread(file_tom,'RAINNC')   )  ) + rot90(  double(   ncread(file_tom,'I_RAINNC')   )  )*50;	% 50 is the bucket value in WRF, please be consistent with namelist
            RAINNC_aft = RAINNC_tom - RAINNC_tod;
            
            W_Pre_bef(:,:,n,y) = RAINNC_bef;
            W_Pre_aft(:,:,n,y) = RAINNC_aft;

        end
		
        disp( [num2str(time_year(y)), '-', num2str(m,'%02d')] )

    end

end

clear Pre_tod Pre_tom Pre_bef Pre_aft

toc 






